Quantum thermometry using the ac Stark shift within the Rabi model 
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A quantum two level system coupled to a harmonic oscillator represents a ubiquitous physical system. New 
experiments in circuit QED and nano-electromechanical systems (NEMS) achieve unprecedented coupling 
strength at large detuning between qubit and oscillator, thus requiring a theoretical treatment beyond the Jaynes 
Cummings model. Here we present a new method for describing the qubit dynamics in this regime, based on an 
oscillator correlation function expansion of a non-Markovian master equation in the polaron frame. Our tech- 
nique yields a new numerical method as well as a succinct approximate expression for the qubit dynamics. We 
obtain a new expression for the ac Stark shift and show that this enables practical and precise qubit thermometry 
of an oscillator. 



The qubit-oscillator model has gone by many names in 
many fields, owing its tenacity to the breadth of its applica- 
bility: It is the simplest non-trivial model of the interaction 
between light and matter. At its inception it was used to de- 
scribe the interaction of an atom with a magnetic field [l|], 
and referred to thereafter as the Rabi model. In the subse- 
quent decades it has been extensively studied in quantum op- 
tics fl2[] and cavity QED [3fl. Physical chemists have used 
a 'vibration-dimer' model to study the spectra of molecules 
|4[] . Applying the rotating wave approximation (RWA) to the 
Rabi model yields the Jaynes Cummings model (JCM) ]5|], 
which is valid when the detuning between the qubit transi- 
tion frequency ft and the resonator frequency u> is negligible 
(n rs w) and the coupling between the qubit and oscillator is 
weak (g < u>) |6J|. This is an excellent approximation in the 
case of cavity QED where typical coupling strengths are of 
order g/u) ps 10~ 6 . The JCM can be extended to incorporate 
tunnelling, and has provided an adequate description of exper- 
iments for decades, but a new era of experiments are pushing 
beyond its boundaries in terms of both coupling and detun- 
ing. Circuit QED experiments couple superconducting qubits 
to LC and waveguide resonators, allowing coupling strengths 
up to g/uj ps 10 -1 , recently enabling demonstrations of the 
breakdown of the JCM [7|, 18|]. Superconducting qubits cou- 
pled to nanomechanical resonators (NR) generally have more 
modest coupling strengths |9l LLOfl, but combined with large 
detuning they could also operate outside the validity of the 
JCM El. 

The Hamiltonian for the Rabi model can be decomposed 
into three parts: 

H = H Q + H + H I . (1) 

The qubit, atom or two level system is described by: 

e A 

where a z and a x are the Pauli spin operators. They describe a 
two level system with an energy splitting e and a spontaneous 
tunnelling between the states at a rate A. In isolation such 



a system would undergo Rabi oscillations with a frequency 



= V^ 



A 2 . The Hamiltonian of the oscillator is: 

Ho = uia'a, 



(3) 



where u is the frequency of the oscillator and a) and a are 
its creation and annihilation operators respectively. Note we 
have neglected the zero point energy. The Hamiltonian for the 
interaction between the two is: 



Hi = g(a + a')a z , 



(4) 



where g is the coupling strength between the qubit and oscil- 
lator. 

Recent experimental progress has sparked a renewed theo- 
retical interest in extending solutions of (Q3 beyond the RWA. 
For instance, a change of basis prior to applying the RWA 
leads to a generalised RWA that should be valid outside the 
very weak coupling limit 11211 . However, this is limited to the 
case of e = 0. As an alternative approach, Van Vleck pertur- 
bation theory fl 1 311 has been used to investigate the dynamics 
in the ultra strong (g/u) > 1) coupling regime Ill4lll5ll . This 
approach contains the splitting and tunnelling elements, but it 
is perturbative in the latter and fails to recover the JCM in the 
weak coupling limit. This approach is therefore more appli- 
cable to circuit QED, rather than the more modest couplings 
achieved in Cooper pair box (CPB) coupled to NR systems. 

An analytic expression for the eigenspectrum of the full 
Rabi model was very recently found by Braak 01611 . a surpris- 
ing and significant result for such a long standing problem. 
However, it is too early to tell how much this solution can 
reveal about the physical properties and dynamics of the sys- 
tem. In addition to solving the model Braak proved that it is 
non-integrable i.e. the time-dependence of important proper- 
ties cannot be found in closed form. There is therefore still 
a need for approximate results governing areas of particular 
experimental interest. 

In order to simplify the expression and extend the validity 
of the approximations that we will subsequently describe into 
the strong coupling regime, we first perform a 'polaron' trans- 
formation IU7U18I1 . This unitary Hamiltonian transformation 



(H' = e s He~ s ) is equivalent to dressing qubit excitations 
with the vibrational modes to form quasi-particles called po- 
larons. With s = a/2(a j( — a)a z and a/2 = g/cu we obtain 

H' = 1<t z + «a t a + ^(D(a)|0)(l| + D(-o)|l)<0|), (5) 

where -D(£) = cxp(£a^ — £*a) is the displacement operator. 
We have neglected a term proportional to the identity g 2 /uj 1, 
which does not influence the dynamics. The first two terms 
involve the qubit and oscillator individually and so can be re- 
moved by going to the interaction picture. We insert the re- 
sulting Hamiltonian into the von Neumann equation and then 
derive equations of motion for the qubit (see Appendix and 
Refs. CI 3): 



^PooW = -»t(pio(«) - Poi(t))), (6) 



&Pn(*) 



.A. 
2 

*2-(Pio(*)-Poi(*))). (7) 



where poo(t) and p\\ (t) are the time dependant population el- 
ements of the qubit's reduced density matrix. The coherences 
are given by: 

P oi(t) = if Jo dt'e-^ [p 00 (f)C*(r) - pii{f)C(r)] ,(8) 
Pio(t) = -t# /J df e*" [a>o(*')C(t) - pu (*')<?* (r)] ,(9) 

where r = i — i' and C(t) and C* (r) are the correlation func- 
tion of the oscillator and its complex conjugate respectively. 
In deriving these equations, we have employed the Born ap- 
proximation, i.e. we have assumed that the vibrational mode 
and the qubit states can be factored at all times. 

Physically, this corresponds to an oscillator that thermalises 
on a timescale faster than that characteristic of the qubit dy- 
namics. 

The equations of motion take the form of a system of 
integro-differential equations involving the bosonic correla- 
tion function and its complex conjugate. Laplace transform- 
ing the equations of motion yields a set of simultaneous equa- 
tions that can be solved algebraically: 



^oo(s) 

Rw(s) 



s* + s 

A " 



spo + (±) 2 [C' + + C>L] 

(10 2 [cl + ci + c + + c;] 



(C'_ + Cl)Roo(s) - -C! 



(10) 



(11) 



where s is our Laplace space variable, i?oo( s ) an d Rio{s) are 
the Laplace transforms of poo(i) and pw(t), Po is the initial 
population of the ground state and C'± = C'(s ±ie) and C± 
are the Laplace transforms of the correlation function and its 
conjugate respectively. It is sufficient to solve these two equa- 
tions alone because from their solutions the behaviour of the 
other density matrix elements can be trivially derived. 

To obtain expressions for the dynamics of Eqns ([Tot and 
( fTTT i in the time domain we need to find the Laplace trans- 
form of the bosonic correlation and its conjugate, solve and 



then take the inverse Laplace transform of the equations. The 
correlation function is defined as 



C(t) = (D t (a)Dl(a))=Tr B \p B D t (a)D\,(a) 
which evaluates to (see Appendix): 

fU T \ _ p— |a| 2 ((l-cos (ui)) coth ^i+jsin (ur)) 



(12) 



(13) 



Unfortunately, it is not straightforward to Laplace transform 
this expression directly, so we employ the Jacobi- Anger series 
expansion: 



E *«(*)« 



i n() 



(14) 



where z is an arbitrary complex number and /„ (z) is the mod- 
ified Bessel function of order n and argument z. By exploiting 
an angle addition identity we can rewrite (fT3~t as: 



C(r) 



— |ck| coth (■^~) „z cos (wr+cc) 



where x = i0w/2 and z = 2\af^/N(N + 1) 
Eqn (O this gives: 



C(r) 



\a\ 2 (2N+ 



1} E J »(*) e 



in(oJT-\-x) . 



(15) 

Using 

(16) 



where N = (e^ u — 1) _1 is the average oscillator occupation 
number. In this form the correlation function can be Laplace 
transformed trivially. The Bessel function weighting of the se- 
ries means it converges very rapidly, in fact only a few terms 
of the series need to be retained to accurately capture the dy- 
namics. For experimentally relevant parameters (i.e. low tem- 
peratures and moderate to strong coupling) only a single term 
dominates. The corresponds to the regime where: 



2\a\ 2 y/N(N + 1)«1. 



(17) 



Physically this corresponds to retaining only interactions 
that conserve the total phonon number in the oscillator, thus 
complementing the underlying Born Approximation, which 
assumes the oscillator remains in thermal equilibrium. In- 
cluding only the dominant zeroth term in the series allows the 
equations ( fT0T > and (fTTT > to be inverse Laplace transformed: 
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+ i&.sm(tri) - e) 
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n=JA 2 e- b I (z) + e 2 1 



(19) 
(20) 



where b — |a| 2 (2iV + 1). From Eqn (|20] | we can see that 
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FIG. 1: Comparison of the single term approximation (red) and a 
numerically exact approach (blue) for different coupling strengths. 
Uncoupled Rabi oscillations are also shown as a reference (green). 
Left: the population poo (t) in the time-domain. Right: the same data 
in the frequency domain. The full numerical solution was Fourier 
transformed using Matlab's FFT algorithm. Other parameters are 
uj = 1 GHz, g = 0.1 GHz, e = A = 100 MHz and T = 10 mK. 




.. L 


L 



0.1 0.2 0.4 0.6 O.f 

CO 



100 







i 0.8 


o 






II 




kp.6 


V 










0\4 








u 




o\ 



1 -3-2-10123 

n 



FIG. 2: Main panel: comparison of dynamics calculated from trun- 
cating l |16| > at JVmax = ±10 (red) and a numerically exact approach 
(blue). Lower left: Fourier transform of the dynamics. Lower right: 
the numerical weight of the n th term in the series expansion of ( 1161 , 
showing there are still only two dominant frequencies at n = 
and n — — 1. Parameters: uj = 0.5, g = 0.1, e = 0, A = 0.5, 
T = 1 mK, h = 1 and k b = 1. 



the presence of the oscillator alters the tunnelling rate by a 
factor Io(z)e~ b , essentially corresponding to a temperature 
dependent ac Stark shift on the qubit due to the presence of the 
single oscillator mode lull 12211 . In contrast to previous work 
our expression is not confined to the weak coupling or large 
detunning limit, but still takes a surprisingly simple closed 
form. 

FigureQ]shows a comparison of the dynamics predicted us- 
ing these expressions and a numerically exact approach. The 
latter are obtained by imposing a truncation of the oscilla- 
tor Hilbert space at a point where the dynamics have con- 
verged and any higher modes have an extremely low occu- 
pation probability. Our zeroth order approximation proves to 
be unexpectedly powerful, giving accurate dynamics well into 
the strong coupling regime (g/u = 0.25) and even beyond 
this it still captures the dominant oscillatory behaviour, see 
Figure Q] Stronger coupling increases the numerical weight 
of higher frequency terms in the series, causing a modulation 
of the dynamics. The approximation starts to break down at 
(g/oj = 0.5). The equations ( U~8T l and ( [T9l ) are obviously un- 
able to capture the higher frequency modulations to the dy- 
namics or any potential long time phenomena like collapse 
and revival, but these are unlikely to be resolvable in experi- 
ments in any case. Nonetheless, it is worth pointing out that 
even in this strong coupling case the base frequency of the 
qubit dynamics is still adequately captured by our single term 
approximation. 

Our methodology can be used to predict dynamics of 
nanomechanical resonators connected to either quantum dots 
or superconducting qubits. The criterion for the single term 



approximation to be valid is readily met by current experi- 



et by ci 



ments such as those presented in Refs. ]9|, |10|] and their pa- 
rameters yield near perfect agreement between numerical and 
analytic results. Most experiments operate in a regime where 
the qubit dynamics are not greatly perturbed by the pres- 
ence of the oscillator, which has a much lower frequency 
(e w A w 10 GHz, u = 1 GHz). In Figure [Q we chose 
e « A w 100 MHz, because this better demonstrates the ef- 
fect of the oscillator on the qubit. These parameters can be 
achieved experimentally using the same qubit desi gn b ut with 
an oscillating voltage applied to the CPB bias gate 11 111 . How- 
ever, we stress the accuracy of our method is not restricted to 
this regime. 

Including extra terms in the series expansion (Q~6j makes the 
time dependence of the qubit dynamics analytically unwieldy, 
because the rational function form of the series leads to a com- 
plex interdepence of the positions of the poles in (101 . How- 
ever, if the values of the parameters are known the series can 
truncated at (±7Vmax) to give an efficient numerical method 
to obtain more accurate dynamics, extending the applicability 
of our approach beyond the regime described by (ITTi . This 
is demonstrated in Fig. [2] where the dynamics are clearly 
dominated by two frequencies - an effect that could obviously 
never be captured by a single term approximation. There is a 
qualitative agreement between the many terms expansion and 
full numerical solution, particularly at short times. We would 
not expect a perfect agreement in this case because the sim- 
ulations are of the dynamics in the large tunnelling regime 
(A = 0.5), and the polaron transform makes the master equa- 
tion perturbative in this parameter. Nonetheless, the rapid con- 
vergence of the series is shown in Fig. [2] Nmax = 5 — 10 is 
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FIG. 3: Demonstration of qubit thermometry: Ti n is the temperature 
supplied to the numerical simulation of the system and T ou t is the 
temperature that would be predicted by fitting oscillations with fre- 
quency J20t to it. The blue line is the data and red line shows the ef- 
fect of a 10kHz error in the frequency measurement; the grey dashed 
line serves as a guide to the eye. The lower inset shows the variation 
of the qubit frequency O with temperature. The upper inset shows 
the dependence of the absolute error in the prediction against the sig- 
nal length (see text). Other parameters are: u) — 1 GHz, g = 0.01 
GHz, e = 0, A = 100 MHz 



sufficient to calculate poo(t) and pio(t) with an accuracy only 
limited by the underlying Born Approximation. The asymme- 
try of the amplitudes of the terms in the series expansion of 
( [Tol l is due to the exponential functions in the series. 

We now discuss the application of our model system to the 
measurement of the temperature of an oscillator by observ- 
ing the coupled qubit. We picture a situation in which the 
qubit oscillation frequency O is the measured quantity. The 
tunnelling, coupling strength and energy splitting are usually 
within the control of the experimentalist (or are at least known 
constants), and this yields the possibility of using the mea- 
sured il to estimate the temperature. A related idea was re- 
cently used in the calibration of a seminal resonator experi- 
ment [| 1 Oil to verify that the oscillator was in its ground state 
(a critical part of the work). In that case, the authors used a 
comparison of numerical results for different occupation num- 
bers N with the measured population in the excited state of the 
qubit after a certain interaction time. A theoretical study of the 
same approach was preformed in H23I1 , where the system was 
described by the JCM without a tunnelling term. In contrast, 
we here propose that simple analytic expressions for the qubit 
dynamics that are valid beyond the weak coupling regime, like 
our own, can be used to directly measure the temperature and 
hence N of the oscillator, simply by observing the effective 
qubit Rabi frequency Q. 

Figure [3] demonstrates this idea, showing that by mea- 



suring D, and fitting it to our expression (120> . we can ob- 
tain submilli-Kelvin precision in the experimentally relevant 
regime of 20-55 mK. At low temperatures the single term fre- 
quency plateaus, causing the accuracy to break down. In the 
higher temperature limit, we also see a deviation from the di- 
agonal, this is to be expected as we leave the regime of valid- 
ity described by $1% . Naturally accuracy in this region could 
be improved by retaining higher order terms in (Q~6}, but this 
would become a more numeric than analytic approach. The 
upper inset shows the dependence of the accuracy of the pre- 
diction on the number of points (at a separation of Ins) sam- 
pled from the dynamics. The accuracy increases initially as 
more points improve the fitted value of Q, however after a 
certain length the accuracy is diminished by long term enve- 
lope effects in the dynamics not captured by the single term 
approximation. We note that the corresponding analysis in 
the frequency domain would not be equally affected by the 
long time envelope, however a large number of points in the 
FFT is then required in order to obtain the desired accuracy. 
The lower inset of Figure [3] shows the direct dependence of 
O on the temperature. The temperature range with steepest 
gradient and hence greatest frequency dependence on temper- 
ature varies with the coupling strength; thus the device could 
be specifically designed to have a maximal sensitivity in the 
temperature range of the most interest. 

In summary, we have developed and explored a new ap- 
proach to the qubit-oscillator model, finding succinct expres- 
sions for the qubit dynamics. In contrast to previous theo- 
retical approaches, our expressions are valid in the stronger 
coupling regime that is rapidly gaining experimental rele- 
vance. We have further proposed an application of our model 
enabling precise temperature measurements of the oscillator 
mode. This could be used either as part of the calibration of 
an oscillator experiment or as a tuned, standalone device. 
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Appendix: derivation of the equations of motion and the bosonic correlation function 



In this Appendix we give the explicit derivation of the equations of motion (|6]|9j, the Laplace transformed equations of motion 
( I10I1U . and the bosonic correlation function ([T2l i. We note that parts of these derivations can be found in similar form in the 
literature (cf. Refs. IU8W20I1 ). but we here give an alternate and full account in consistent notation for the benefit of the reader. 



Equations of Motion 



First we move into the interaction picture: 



-lit 7-)t 



Poo(i) = Poo, Pn(t) = pn, Poi(t) = Ane I£ D t , p w (t) = p w e le Dj, 
where D t and Dj are the time dependent versions of the displacement operators introduced by the polaron transform. 

Hi (t) = - (poi(t) + Pw (t)) , 



(21) 



(22) 



Starting from the Von Neumann equation: 

&p(t) = -*[Hi(t),m], (23) 

p(t) = p -iJ*dt'[H I (t'),p(t% (24) 

To study the dynamics we need the time dependent expectation values of the density matrix elements. These are given by: 



Substituting in ( f24l > 



(0) t =Tr[p(t)0}=Ti[p(t)Ot], 



(0) t - (O)a = -i / dt"&[[irj(f),/5(t , )]Ot] 1 

Jo 



Exploiting the cyclic property of traces: 



(25) 



(26) 



(0) t -{O) = -i dt'TiWftOuHrf)]], 
Jo 



(27) 



Substituting O for the relevant operator eg. poo(t), evaluating the commutator, and tracing over the qubit degrees of freedom 
yields: 



(poo(*)> - (poo(O)} = -*# /„* dt'((p 10 (t')} - (poi («')», (28) 

(pn(*))-(pn(0))= if f>'((p 10 (£')>- W (*')>), (29) 

(poi(*)} - (Poi(0)) = -if ^^-^((poo^A^) - (pii(f)4A», ( 3 °) 

(pio(t)) - (pio(O)) = i#/o ^-"('-^((^(OA^I) - (pii(*0A t A'»> (3D 

At this point we make the Born approximation (assuming the density matrix of system and bath are factorable) 

(p o(t')D t (a)D t ,(a))t> « ( Poo (t'))(D t (a)D t ,.(a)) (32) 
The bosonic correlation function is defined as C(t — t'): 

C(t - t') = {D t {a)D\,{a)) = Tr B [p B D t {a)D\,{a)\. (33) 



Where the subscript B represents the bosonic degrees of freedom. We substitute this into (128b and by assuming there is no initial 
coherence in the system we obtain (|6]-|9j. 

Correlation Function 

The bosonic correlation function (fT2l for an oscillator with a single mode in a thermal state is defined as: 

C(t t') = Tr B [p B D t (a)Dl(a)}, p B = -^tf^l = )_ C xp(-/Wa). (34) 

This can be evaluated in different ways, one of which is presented below. Starting from the time dependence of the displacement 
operator in the interaction picture: 

D t (£) = e iHot D(C)e~ iHot = e iuia Klt D(£)e~ i " a >at , (35) 

or, alternatively, through the time dependence of creation and annihilation operators: 

D t (£) = e^^-Cae-^ = D £ e *uty (36) 



In order to perform the trace Trs in the number state basis, we need to know the action of e^ a a and -D(£) on a number state 
\n). The first simply evaluates to e^ n and the latter gives the so-called displaced number state |£, n). 
The displaced number state can be expanded in the number state basis 

oo 

\tn) = Y, C nm \m), C nm - (m|D(0|n), (37) 

m=0 



with (see, e.g., Oliviera et al 12411 or M. Crisp I125I1 ) 



c nm = \—,e-^ 2 r~ n u;r n m 2 ), (38) 

V ml 

where £™ -n (|£| ) is an associated Laguerre polynomial. This is only valid for m > n, but for m < n the displacement operator, 
or rather its hermitian conjugate, can be made to act on (m\ instead of on \n). 

Single series and analytical result 

We use Eq. ( I36J I for the displacement operator and the property D(x)D(y) = cxp[(xy* — yx*)/2]D{x + y) to evaluate Eq. 
(1341 . This leads to a series of the following form 

C(t - t') = _ e -M a [ 1 - e ~ Mt ~ t ' ) ] Y, e- 0un L n [2\a\ 2 (l - cos[oj(t - t')})}. (39) 



7 

By virtue of the property J2n=o L n (y)z n = (1 - z)' 1 exp[yz/(z - 1)] and with N = (e^ - l)" 1 and Z = (1 - e-^y 1 we 
finally arrive at the analytical result: 

C(t - if) = cxpHH 2 sin[w(i - if)]} exp[-2|a| 2 (l - cos[w(i - t')]){N + 1/2)]. (40) 

Note that this expression agrees with Mahan's result for a single mode (Ref 11811 . section 4.3) C(t) = 
e -|a| 2 ((i-cosa;t)coth(^)+i 8 in^ Mahan derives this in a similar fashion but without using Eqs. 071 ED- Instead, he uses 
the 'Feynman disentanglement of operators' to arrive at an equivalent infinite series of Laguerre polynomials. 



